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Abstract 

This is the first of two papers about the structure of Kauffman networks. In this 
paper we define the relevant elements of random networks of automata, following previ- 
ous work by Flyvbjerg Q and Flyvbjerg and Kjaer and we study numerically their 
probability distribution in the chaotic phase and on the critical line of the model. A 
simple approximate argument predicts that their number scales as ^fN on the critical 
line, while it is linear with N in the chaotic phase and independent on system size 
in the frozen phase. This argument is confirmed by numerical results. The study of 
the relevant elements gives useful information about the properties of the attractors 
in critical networks, where the pictures coming from either approximate computation 
methods or from simulations are not very clear. 

1 Introduction 

Kauffman networks are disordered dynamical systems proposed by Kauffman in 1969 as a 
model for genetic regulatory systems |1|. They attracted the interest of physicists in the 
80's P, 1^, H, H, due to their analogy with the disordered systems studied in statistical 
mechanics, such as the mean field Spin Glass [0. A dynamical phase transition was found 
and studied in the framework of mean field theory. 

In this and in the next paper [^] we deal with some structural properties of the networks 
that determine their attractors. In the present paper we introduce the relevant elements, 
a notion that was suggested by Flyvbjerg p and Flyvbjerg and Kjaer [^, and we study 
their probability distribution. In the next one we describe how the relevant elements are 
subdivided into asymptotically non communicating, independent modules. The modular 
organization of random boolean networks was already suggested by Kauffman IQ , and it was 
used by Flyvbjerg and Kjaer to study analytically the attractors in = 1 networks. We 
shall show that it is possible to describe the phase transition in random boolean networks 



in terms of the scaling of the number of relevant elements with system size, or in terms of 
a percolation transition in the set of the relevant elements. The interest of this approach 
is that some consequences about the statistical properties of the attractors can be directly 
drawn. 

In IP we computed the properties of the attractors in the framework of the annealed 
approximation, introduced by Derrida and Pomeau [0, but we observed that the results of 
this approximation are reliable only when the system is chaotic enough, becoming exact for a 
random map. The study of the relevant elements is complementary to this approach, and we 
sketch the lines of a new approximation scheme that works better in the frozen phase and on 
the critical line. This region in parameter space is the most interesting one, since, according 
to Kauffman, it reproduces some features of real cells, and is also the less understood, since 
neither approximate computations nor simulations []TU[ give a precise picture of the properties 
of the attractors for systems of large size. 

In next section we define the model, discussing some old results together with open prob- 
lems. In section 3 we define the relevant elements and in section 4 we give an approximate 
argument predicting the scaling of their number with system size in the different phases 
of the model. In the following section we present our numerical results, starting from the 
magnetization and the stable elements (section 5.1) and then discussing the distribution of 
the relevant elements and its connection with the properties of the attractors, respectively 
in the chaotic phase (section 5.2) and on the critical line (section 5.3). The discussion of the 
results is postponed to our following paper [Q, concerning the modular organization of the 
relevant elements on the critical line. 



2 Definition of the model and previous works 

Kauffman model is defined as follows. We consider a set of elements Q = {1,---A^} 
and we associate to each of them a binary variable, cTj G {0, G Q. In the biological 
interpretation proposed by Kauffman each element of the network represents one gene and 
the binary variable cTj represents its state of activation. 

Each element is under the control of K elements, in the sense that its state at time 
t + 1 is determined by the states at time t of the K control genes, ■ ■ ■ ixii) and by 

a response function of K binary variables, /j(o"i, ■ ■ ■ o"/<) G {0,1}, that specifies how the 
element i responds to the signals coming from its control variables. The control elements are 
chosen in f2 with uniform probability. The response functions are also extracted at random, 
and it's believed that the properties of the model do not depend on the details of their 
distribution |^, |T0[. The rule most generally used in teh literature is the following: for 
each of the possible inputs G {0, 1}^ we extract independently the value of and we call 
p the probability that is equal to 0. 

The dynamics of the system obey the equation 

ai{t + 1) = /i ((Tji(i), ■ ■ ■ (Tj^i^i)) . (1) 

This evolution law is deterministic, but the system is disordered because the control rules 
(elements and functions) are chosen at random from the beginning and kept fixed: thus we 
deal with a statistical ensemble of deterministic dynamical systems, and we are interested 



in the statistical properties of systems of large size. For finite N, every trajectory becomes 
periodic after a long enough transient time, and the configuration space is partitioned into 
the attraction basins of the different periodic orbits. We are interested in the probability 
distributions of the number, the length and the size of the attraction basin of the periodic 
orbits, as well as in that of transient times. In the biological metaphor, given a set of rules (a 
genome) an attractor represents a possible cellular type, its length represents the duration 
of the cellular cycle, and the number of attractors represents the number of cells that can 
be formed with a given genome. 

It was observed already in the first simulations that two dynamical regimes are present, 
and that the line separating them has properties reminiscent of those of real cells [0]. In the 
so-called chaotic phase (large connectivity, p close to 1/2) the average length of the cycles 
increases exponentially with system size. The limit case of the chaotic phase, K — >■ oo, was 
already known as Random Map in the mathematical literature, and was studied in detail 
by Derrida and Flyvbjerg |jTl|, who pointed out interesting analogies between this system 
and the mean field Spin Glass |^ concerning the distribution of the weights of the attraction 
basins. In the frozen phase, on the other hand, the typical length of the cycles does not 
increase with N. The limit case of this phase, K = 1, was analytically studied to some 
extent by Flyvbjerg and Kjaer 0, who introduced in that context the concept of relevant 
elements (though without using this name). 

The first description of this dynamical phase transition in terms of an order parameter 
was given by Derrida and Pomeau 0. They studied the evolution of the Hamming dis- 
tance between configurations in the Kauffman networks approximating it with a Markovian 
stochastic process. Such approximation (the so-called annealed approximation) was then 
shown to be exact in the infinite size limit, concerning the average value of the distance 0. 
Below a critical line in parameter space the average distance goes to zero in the infinite size 
limit {frozen phase) and above it the distance goes to a finite value {chaotic phase). The 
position of the phase transition depends only on the parameter p, representing the proba- 
bility that the responses to two different signals are different^], and is given by the equation 
p,{K) = 1/K. 

The properties of the attractors can be easily computed from the knowledge of the whole 
stationary distribution of the distance, and this can also be obtained within the annealed 
approximation j^, but the validity of this approximation in this more general case is not 
guaranteed. Comparison with simulations shows that the agreement is satisfactory in the 
chaotic phase, while the approximation fails on the critical line. In the chaotic phase it 
is possible to compute the value of the exponent of the typical length of a cycle, r cx 
exp (a(i^, p)A^), in good agreement with numerical results, but the distribution of cycle 
lengths is much broader than it is expected. The annealed approximation predicts also that 
the distribution of the weights of the attraction basins is universal in the whole chaotic 
phase, and equal to the one obtained by Derrida and Flyvbjerg in the case of the Random 
Map The corrections to this prediction appear small, if any, even for K = 3. Finally, 
the number of different cycles in a network is expected to be linear in N, but it is very hard 
to test numerically this prediction. 

The annealed approximation makes also predictions about the critical line of the model 



^In terms of p one has p — 2p{l — p), so its value is comprised between zero and 1/2, but for K = 1 p 
can be taken as an independent parameter in [0,1] 



0. It predicts that the properties of the attractors are universal on the critical line p = 1/K 
(with the exceptions of the points K = 1, p = 1 and K = oo, p = 0, which are not transition 
points). In particular, the typical length of the cycles should increase as all along the 
critical line. Numerical results are not clear under this respect [0: it seems that the rescaled 
cycle length / = L/\/N has a limit distribution if / is small (roughly, smaller than 2) but 
for larger values the distribution becomes broader and broader as increases ||T^ , p!0| , so 
that it is possible to define an effective length scale increasing much faster with system size 
(as a stretched exponential). The distribution of the number of cycles has exactly the same 
characteristics. These results cast doubts on the validity of the biological analogy proposed 
by Kauffman, that relies very much on the fact that in critical networks the typical number 
of cycles scales as V^, reminiscent of the fact that the number of cell types of multicellular 
organisms very far apart in the filogenetic tree scales as the square root of the number of the 
genes, and that in critical networks the typical length of the cycles increases as a power law 
of system size, also consistently with the behavior of cell cycles time. Thus it is interesting 
to understand how these distributions look like in the limit of very large systems. 

Another reason of interest of the present approach is that it allows to understand the 
limits of the annealed approximation. In our interpretation the annealed approximation 
is valid as far as the system loses memory of the details of its evolution. This, of course, 
does not happen if in a realization of a random network some structural properties that 
are able to influence its asymptotic dynamics emerge. Thus the approach presented here is 
complementary to the one used in 0. 

3 Definition of the relevant elements 

Let us start recalling the definition of the stable elements 0. These are elements that evolve 
to a constant state, independent of the initial configuration. Flyvbjerg defined them and 
computed their fraction s = S/N using the annealed approximation, which becomes exact 
in the infinite size limit. We now recall briefly, for future convenience, the main steps of this 
calculation. 

Let us suppose that an element is controlled hjK — i stable elements and i stable ones. 
Then it will be stable if the control function does not depend on the unstable arguments 
when the stable arguments assume their fixed values. Otherwise it will be unstable. When 
all the i unstable elements are different (this can always be taken to be the case if K is 
finite and grows), the probability Pi to choose a constant control function of i binary 
variables is given by Pi = p"'* + (1 — p)"S with ui = 2\ In the framework of the annealed 
approximation, extracting at random connections and response functions at each time step, 
we get the following equation for the fraction of variables that are stable at time t: 

s{t + 1) = 7 {s{t)) = E f f ) s{tf-^ (1 - s{t)y P,. (2) 

This equation can be shown to be exact in the infinite size limit. The fixed point of 
this map (which can be interpreted as a self-consistency equation for the fraction of stable 
variables) has only the trivial solution s = 1 in the frozen phase, in other words all the 
elements are stable except eventually a number increasing less than linearly with A^. In the 



chaotic phase this solution becomes unstable and another solution less than 1 appears. This 
happens when K{1 — Pi) = 1. Since 1 — Pi = p (it is just the probability that the response 
to two difTerent signals are different) this condition is equivalent to the condition obtained 
from the study of the Hamming distance. 

The existence of the stable variables is due to the finite connectivity of the network (s* 
goes to zero very fast when K increases) . These variables do not take part in the asymptotic 
dynamics. Among the remaining unstable variables, some are irrelevant for the dynamics, 
either because they do not send signals to any other variable, or because they send signals, 
but the response functions are independent of this signal when the stable variables have 
attained their fixed values. The remaining variables, that are unstable and control some 
unstable variable, are what we call the relevant variables. They are the only ones that can 
influence the long time behavior of the system. 

To be more clear we now describe the algorithm that we used to identify the relevant 
variables. As a first step, we have to identify the stable variables. These are the variables that 
assume the same constant state in every limit cycle, and identifying them is computationally 
very hard, but very simple in principle. We then eliminate from the system the stable 
variables, reducing the response functions to functions of the unstable variables alone. Some 
of the connections left are still irrelevant, and we have to eliminate them (a connection 
between the elements i and j is irrelevant if the reduced response function fi{(Tj-^{i)-i ■ ■ ■ crjidi)) 
does not depend on the argument for all the configurations of the remaining Ki — 1 control 
variables). At this point we iterate a procedure to eliminate the irrelevant variables. At each 
iteration we eliminate the variables that do not send any signal to anyone of the variables 
that are left, until we remain with a set that cannot be further reduced. This is the set of 
the relevant variables. 

Measuring the number of relevant variables is computationally a very hard task. In order 
to identify the stable variables, in fact, we should find all the cycles in the network, and, to 
be rigorous, we should simulate a number of trajectories of the same order of the number 
of configurations in the system. Of course this is not feasible and we run only 200 (in some 
case 300) randomly chosen trajectories in every network. Thus we overestimate the number 
of stable elements. Nevertheless, the number of stable elements changes very little when we 
simulate more initial conditions and we think that the error that we make is not very large. 
However, for every network we simulate some hundreds of trajectories and every trajectories 
has to be followed until the closing time. This grows exponentially with system size in the 
chaotic phase. On the critical line the typical closing time increase roughly as a power law of 
system size, but the distribution becomes broader and broader and the average closing time 
is more and more dominated by rare events. The average depends thus on the number of 
samples generated and on the cutoff of the closing time, i.e. the maximum time that we are 
disposed to wait to look for a cycle. To reduce the bias determined by the cutoff, we had to 
run simulations lasting a time which increases roughly as a stretched exponential of system 
size on the critical line. Thus it is not possible to simulate systems of more than about one 
hundred elements in the chaotic phase and one thousands of elements on the critical line. 



4 Scaling argument in the frozen phase 



The mean field analysis shows that the fraction of relevant variables vanishes in the frozen 
phase and on the critical line, but does not tell how the number of relevant variables scales 
with N as N grows. In order to clarify this point, we have to go beyond the mean field 
picture. 

In the special case of i^' = 1, belonging to the frozen phase for every p < 1, there are 
detailed analytical results about the distribution of the relevant variables 0. We propose 
here a rough argument that generalizes those results to the whole frozen phase and predicts 
that the typical number of relevant elements scales as \/N on the critical line. Though this 
argument is based on some approximations which we can not control, its results coincide for 
K = 1 with the exact results by Flyvbjerg and Kjaer. 

Let us suppose that we add a new element to a system with N elements, R of which are 
relevant, while S are stable and / = N — R—S are indifferent, i.e. neither stable nor relevant. 
The probability that the new element is relevant can be computed as a function of R and 
S, within some approximations that we are going to discuss in a while. This probability is 
equal to the fraction of relevant elements in the system with + 1 elements, given that the 
relevant elements are R and the stable ones are S in the system with N elements. We can 
then average over R and S in order to get an equation connecting r^r+i = {R) jy^i/ {N + 1) 
to the moments of the distribution of R in the system with N elements. Since in the frozen 
phase and on the critical line r^v vanishes, it will be enough to consider the first two moments 
of the distribution, and the resulting equation can be solved asymptotically in A^. 

The weakness of this approach lies on the assumptions that allow us to express the 
probability that the new element is relevant as a function of R and S, as it will become soon 
clear. We compute now this probability. To this aim, we need two steps: 

1. As a first step, we have to extract the K control elements and the response function 
of the new element. As a consequence, the new element can be stable, unstable or, 
if it receives an input from itself and this input is relevant in the sense discussed 
above, relevant. The evaluation of the stability is perfectly equivalent to the mean 
field argument, but this stability is only temporary because it can be altered by the 
second step described below. Thus we call a new element that is stable (unstable) after 
the first step a temporarily stable (unstable) element. 

2. Then we have to send to the old system the signal of the new element. For each of 
the KN old control connections we have a probability 1/(A^ + 1) that the connection 
is broken and the old control element is substituted by the new element. This step 
perturbs the elements that control the new element and modifies its temporary stability. 
We have no chance to take this into account, unless we use some drastic approximations. 

In the second step, three situations can occur. 

1. If the new element was relevant in the first step, the new step can not modify this 
condition. 

2. If the new element was unstable, it cannot become stable through the feedback of its 
signal. So it will be relevant or indifferent, depending on whether it sends an input to 
at least one relevant element or not. 



3. If the new element was stable, its signal can destabilize some of the elements that 
control it and thus it can become relevant through a feedback mechanism, very hard 
to investigate analytically. 

To compute the probability of case 3, we should know the organization of the network 
in the very detail and not only the number of relevant and stable elements. We propose 
to bypass this difficulty considering a different event: we will consider the new element 
relevant if it receives a signal from a previously relevant element or from itself. This is the 
simplest way to get a closed equation for the average number of relevant elements. In this 
way wc make two errors of opposite sign: on one hand we overestimate the probability that 
a temporarily unstable element becomes relevant, on the other one we underestimate the 
probability that the new element is temporarily unstable and we neglect the probability that 
a temporarily stable element becomes relevant through a feedback loop. 

We think that this method captures at least the qualitative behavior of the number of 
relevant elements. We have then to compare the estimate given by this approximation to 
the simulations, because the approximation is not under control. Wc present this argument 
because its results agree with both the numerical results and with the analytical calculations 
ior K — 1 and because we believe that it is possible to improve this method and to keep the 
approximation under control. 

Since we are interested in the frozen phase, where the fraction of unstable elements 
vanishes in the infinite size limit, we can neglect the eventuality that the new element is 
controlled by more than two elements that were relevant in the old system. The results 
are consistent with this assumption. With these approximations we obtain the following 
equation for the probability that the new element is relevant: 
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where p2 represents the probabihty that a Boolean function of two arguments is not constant 
and in terms of p is given by 
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(4) 



In the frozen phase it is sufficient to consider that the new element receives only one 
signal from the previously relevant elements. So, posing c = Kp, the equation for the new 
fraction of relevant elements, r, is 



{r)N+i^c{r)j^ + —. 



(5) 



The first term represents a new element that receives a relevant signal from one of the 
previously relevant elements, the second term represents a new element that receives its own 
relevant signal. 



Thus the average number of relevant elements is independent on and its asymptotic 
value is 

{R)n = (6) 

This number diverges on the critical line c = 1. In this case, we have to consider also 
the eventuality that the new element receives a signal from two of the previously relevant 
elements. Expanding to the second order in r = R/N, and using the fact that pc = 1/K, we 
get the equation 

««« « {r}„ - (^) (r\ + i (7) 

whence, in the asymptotic regime where the variations of {r)^ are of order r/N, we finally 
get 



-r^U. (8) 

This means that the scale of the number of relevant elements grows, on the chaotic phase, 
as y/N. 

We stress here that these computations are valid because of the finite connectivity of the 
system. If we perform the limit ^ cxd on the above result, we get that the scale of the 
number of relevant elements grow as 1/2^/N. If, instead, we apply the limit K ^ oo prior 
to the limit N —>■ oo we get the trivial critical point p = 0, where all the elements are stable 
after one time step, while for every other p value all the elements are relevant. Thus, the 
two limits do not commute. In fact, the equation (§) for the fraction of stable variables and 
all the computations performed in this section are valid only if we can neglect that the same 
element is chosen more than once to control a given element, i.e. for K <^ N. 

The result (|^) coincides for K = 1 with the analytical computation by Fljrvbjerg and 
Kjaer 0, thus suggesting that the distribution of relevant elements is independent on in 
the whole frozen phase, and depends on the two parameters K and p only through their 
product. This picture agrees with the results of the annealed approximation, which predicts 
that the distribution of the number of different elements in two asymptotic configurations 
is independent on and depends only on the product of the parameters K and p in the 
frozen phase ||^. 

Our simulations confirm that on the critical line the number of relevant elements scales 
as (see figure |]). Also the annealed approximation is consistent with this result, since it 
predicts that the number of elements whose state is different in two asymptotic configurations 
has to be rescaled with ■s/N on the critical line 0. On the other hand the number of unstable 
elements grows much faster with A^ (numerically it is found that it goes as see below) 

but this discrepancy is only apparent, since the asymptotic Hamming distance is related 
more to the number of relevant elements than to this quantity. 

For later convenience (see our next paper) it is also interesting to compute the effec- 
tive connectivity, defined as the average value of the relevant connections between relevant 
elements. Let us compute it by imposing the condition that the network has R relevant 
elements. 

The effective connectivity is equal to the average number of connections between the new 
element and the other relevant elements of the older system, with the condition that the new 
element is relevant. From equation (^ we have, at the leading order in R/N: 



This equation shows that the effective connectivity minus 1 goes to zero as R/N in the 
frozen phase (where R/N cx 1/iV) and on the critical hne (where R/N oc 1/\/N). For a 
fixed system size, the effective connectivity increases linearly with the number of relevant 
elements. 



5 Numerical results 

5.1 Magnetization and stable elements 

As a first step, our algorithm has to identify the stable elements. It does this by measuring 
their magnetization. We thus discuss our numerical results starting from this quantity. 

The magnetization m° of element i on the cycle can be defined as the average activity 
of the element along the cycle: 

< = 7^ E ^^ic)- (10) 

The distribution of this variable, shown in figure Q for K = 3 and = 75, has many 
peaks, corresponding to simple rational values. This perhaps reflects the fact that the 
relevant elements are divided into asymptotically independent modules, so that a cycle can 
be decomposed into several independent shorter cycles. This subject will be further discussed 
in our second paper. 



Our results have to be compared to the analytical work by Derrida and Flyvbjerg |T3 
They defined the magnetization of element i at time t on a given network as the activity of the 
element at time t averaged over many initial configurations and could compute analytically 
its stationary distribution, in the limit N ^ oo, using the annealed approximation, that 
can be shown to be exact for this purpose. The picture they got is different from ours, in 
particular we see peaks much higher than theirs. For instance the peak at | 2m — 1 |= 1, 
which gives information on the size of the stable core of the network, is about 10 times 
larger then expected, and the first moments of the magnetization, that can be computed 
analytically, are larger than the predicted values. Thus we performed other simulations that 
strongly suggest that these discrepancies are finite size effects, and we present an argument 
that explains their origin. 

In order to investigate larger systems we had to change the definition of the magneti- 
zation. The definition ([l0|) is numerically cumbersome, since the measure takes place only 
after that a cycle has been found, and this means, for chaotic systems, that we have to wait 
a time exponentially increasing with A^. Thus we neglected this condition and we measured 
the magnetization of the variable i at time t as the average activity of the variable with 



respect to different initial conditions (this definition coincides with the one used by Der- 
rida and Flyvbjerg). For very large t, when all trajectories have reached a limit cycle, this 
quantity tends to the asymptotic value 

mi = Y.Wc,m'^, (11) 

a 

where is the weight of the basin of cycle Fq, and mf is defined in (^0]). We observed 
that m{t) reaches a stationary value (within some precision) much earlier than the typical 
time at which the trajectories reach their limit cycles. At first sight surprisingly, the time 
after which m{t) reaches its stationary distribution does decrease with system size instead 
of increasing (see figure 

We measured the second and fourth moment of the magnetization in a system with K = 3 
and p = 1/2, and we found a large positive correction to the infinite size values computed 
by Derrida and Flyvbjerg |T^. The values found coincide within the statistical error with 
those obtained from equation (|T0|) for a system of small size for which we did an explicit 
comparison. These values can be fitted to the sum of the infinite size value, that we got from 
Ill3| , plus an exponentially decreasing term. The exponent of the best fit turned out to be the 
same for both the moments that we measured: we found m2{N) 0. 236 + 0. 24-exp(—A^/70), 
and m4{N) ^ 0.128 + 0.26 ■ exp(-iV/70). 

The measure of the magnetization allows to identify the stable elements as the elements 
with J2a Wamf equal either to or to 1. The two definitions of the magnetization gave 
roughly the same number of stable elements in the cases where we could compare the results, 
but with the second method we could consider much larger systems (we recall that the 
difference between the two methods is that in the first case a cycle has been reached while in 
the second one the system is still in some transient configuration). The second method was 
used only to study finite size effects, since it does not allow to identify the relevant elements 
(see below). 

Both the methods overestimate the number of stable elements, since it could happen 
that an element appearing stable in our sample of trajectories (some hundreds) oscillates 
in a cycle that is not reached by any of them. We checked that the results do not change 
qualitatively if we consider a larger number of trajectories. 

The fraction of stable nodes measured in simulations with K = 3 and ranging from 
50 to 200 have been compared to the prediction of the mean field theory by Flyvbjerg. The 
networks with = 50 have a stable core about 10 times larger than the mean field value (in 
this case we measured the magnetization using both the above definitions, while for larger 
systems only equation (|Tl|) was used). The corrections to the mean field value, that is exact 
in the infinite size limit, appear to decay exponentially with a rate identical, within statistical 
errors, to the decay rate of the corrections to the moments of the magnetization: we found 
s(A^) ^ 0.0122 + 0.21 ■ exp(-A^/70). 

For every size of the systems which we simulated the stable core is then much larger than 
it would be in an infinite system. On this ground, we may expect very important finite size 
effects concerning the dynamical properties of the system. 

Summarizing, the distribution of the magnetization for finite systems has the following 



characteristic: 1) The asymptotic value is reached after a time that decreases with system 
size; 2) the corrections to the infinite size values are very large; and 3) these corrections 
decrease exponentially with system size. These apparently strange finite size effects have a 
simple interpretation: they arise as a consequence of the periodic dynamic of the random 
networks. 

The mean field values of the magnetization and of the stable core are computed within 
the annealed approximation without taking into account the fact that the asymptotic dy- 
namic is periodic. As we proposed in |p, the existence of limit cycles must be taken into 
account in the framework of the annealed approximation in this way: if at time t all the 
configurations generated are different {i.e. the trajectory is still open) we treat the quantities 
of interest (distance, magnetization or stable core) as a Markovian stochastic process; if one 
configuration has been found twice (the trajectory is closed) we impose the condition that 
all quantities are periodic. Thus the master equation for the distribution for the number of 
stable variables is, in the framework of the annealed approximation: 



Pr {S{t + 1) = S\ 0{t + 1) I S{t) = S, 0{t)} = j (7(.))^' (1 - 7(^))^"''' (1 - MS, t)) 



Pr {S{t + 1) = S', 0{t + 1) I Sit) = S,0{t)}= (^^' 

Pr {Sit + 1) = S',0it+1) I Sit) = S,0{t)} = 5ss', (12) 

where Sit) is the number of stable elements, s = S/N, 0(t) stands for the condition that 
the trajectory is open at time t (no configuration has been visited twice), 0(t) stands for 
the condition complementary to 0(t) (the trajectory has closed on a previously visited 
configuration) and 7rjv(S', t) is the probability that a trajectory open at time t and with 5* 
stable elements at time t + 1 closes at that time. Finally, 7(5) is given by equation (^. 

We don't know how to compute T^NiS, t), but it is clear that this is an increasing function 
of S for fixed t: the more elements are stable, the more it is likely that the trajectory closes. 
The infinite size value of the stable core is given by equation (|^), that represents the evolution 
of the most probable number of stable variables. It is clear that the corrections to this value 
are positive, and that they go to zero as soon as the closing time becomes much larger than 
the time necessary for the stable core to reach its stationary value in an infinite system 
(where all trajectories are still open). Thus we expect that these corrections vanish as a 
power law of the typical length of the cycles: in the chaotic phase this means that the finite 
size corrections due to this effect vanish exponentially with system size, as we observed 
simulating systems with K = 3 and p = 1/2. Lastly, this argument implies that the time 
after which the distribution of the stable elements becomes stationary is shorter in an infinite 
system than in a small system, where the evolution of Sit) is coupled to the closure of the 
periodic orbits. Thus the correction of the annealed approximation to take into account the 
existence of periodic attractors can account for all the features of the finite size effects that 
we observed. 



5.2 The relevant elements in the chaotic phase 

After having identified the stable elements we detect the relevant elements using the algo- 
rithm described in the second section and we study how this quantity influences the dynam- 
ical properties of the network. The main results are that the average cycle length grows 
almost exponentially with the number of relevant variables in some range of this variable 
and the average weight of the attraction basins has apparently a non monotonic behavior 
versus the number of relevant variables. This qualitative features are the same both in the 
chaotic phase and on the critical line, but the ranges of R in which these things happen are 
quite different in the two cases. We start discussing the situation in the chaotic phase. 

The simulations were done generating at random 20, 000 sample networks and running 
200 trajectories on each of them. The parameters considered in this section are = 3, 
p = 1/2 and system size N ranging from 30 to 60 elements. 

Figure ^ shows the density of the distribution of the fraction r^v of relevant variables, 
tn = R/N. The density relative to the most probable value increases with system size, and 
it appears that tends to be delta-distributed in the infinite size limit, as it is expected 
on the ground of the annealed master equation ([T^). We observe an excess of networks with 
very few relevant elements [i.e. very many stable elements), consistently with the finite size 
effects discussed in last section. This excess seems to disappear in the infinite size limit. 

Then we show the average length of the cycles in networks with R relevant elements 
(figure This quantity increases almost exponentially with R when r = R/N is large, 
while its behavior is different for small r. The crossover takes place at about r = 0.5. Thus 
the number of relevant elements turns out to have a very important influence on the typical 
length of the cycles 

We have also measured the conditional distribution of the length of the cycles in networks 
with R relevant elements. When R is close to the distribution decays as a stretched expo- 
nential with an exponent smaller than one, very close to the one found in the unconditioned 
distribution. Thus the deviation of the unconditioned distribution from the prediction of the 
annealed approximation, that predicts a much narrower distribution, is not a consequence 
of the existence of the relevant elements. 

The other quantity that we measured is the average weight of the attraction basins, Y2, 
defined by the equation 

Y2=j:wi (13) 

a 

where Wa is the attraction basin of cycle a. We used the method proposed by Derrida 
and Flyvbjerg [Q, that is based on the fact that I2 is equal to the probability that two 
trajectories chosen at random end up on the same attractor. 

From our data (not shown) it appears that Y2 has a non monotonic behavior as a function 
of r: for very small r it decreases from the value 1, corresponding to r = 0, reaches a minimum 
and rapidly increases. At large r, Y2 does not seem to be correlated with r (at least within 
the statistical error, that is rather large). We will see in the next paper that the decreasing 
behavior at small r can be interpreted as an effect of the modular organization of Kauffman 
networks. 



5.3 Relevant elements on the critical line 



We simulated systems with K = 4 and the critical value p = 1/4. Systems size ranges from 
120 to 1600. Concerning the statistical properties of the attractors, these networks have a 
behavior very similar to that of the more studied K = 2, p = 1/2 networks ||10|| . 

In these networks, the number of relevant elements appears to scale as viV, in agreement 
with the argument presented in section 4. The number of unstable elements, on the other 
hand, appears to scale as N^^^. This implies that the probability to extract at random an 
element which is relevant, scaling as N~^/'^, is approximately proportional to the square of 
the probability to extract at random an element which is relevant (A^~^/^). 

These scaling laws can be observed both looking at the average quantities and looking 
at the whole distribution. The average number of unstable variables is found to follow the 
power law U oc A^"^, with a = 0.74 ±0.01. We then define the rescaled variable Xu = U/N^^'^, 
and we compare its probability density for various system sizes. As it can be seen in figure 
^ the different curves superimpose within the statistical errors. This suggests that Xu has a 
well defined probability density in the infinite size limit, although our data are rather noisy 
to state this point without doubts. We can distinguish in the distribution three different 
ranges with different characteristics: at vanishingly small values of Xu (ranging from U = 
up to = 4) the density decreases very fast. At intermediate values, roughly up to = 1, it 
looks to decrease approximately as a power law with a small exponent (the best fit exponents 
that we found range from 0.25 to 0.40, showing some tendency to increase with system size). 
Asymptotically, for large x^, the best fit is a stretched exponential, f{x) ~ exp (^—Cx^^, 
with an exponent compatible with (3 = 1.7 ± 1 for all the systems that with studied with A^ 
larger than 240. 

The number of relevant variables was studied in a similar way. Its average value increases 
as a power law of A^, {R) oc R^, with b = 0.52 ± 0.02. The rescaled variable Xr = R/\/N 
looks to have a well defined distribution in the infinite size limit, as it is shown in figure |^, 
where the probability density of x^ is plotted for system sizes ranging from 120 to 1600. For 
large x the density of the distribution is well fitted by a stretched exponential, exp 
with the exponent (3 compatible with the value 13 = 0.56 ± 0.02 for system size larger than 
240. 

The average length of the cycles increases exponentially as a function of the number of 
relevant elements, for r large, and more slowly for r small, just as it happens in the chaotic 
phase. Figure |^ shows on a logarithmic scale the behavior of the average length of the cycles 
as a function of the rescaled number of relevant elements, Xr = R/VN, for different system 
sizes at the critical point K = 4, p = 1/4. 

The average weight of the attraction basins, Y2, has a non monotonic behavior as a 
function of the number of relevant elements, as it happens in the chaotic phase. The value 
of 1^2 (-R) is one for R = 0, then decreases to a minimum value and increases very slowly, as 
it is shown in figure where Y2 is plotted against Xr = R/\fN, for 7^^ = 4,^=1/4 and 
different system sizes. 

Nevertheless, there are two important differences with respect to the chaotic phase: first, 
the range where 1^2 (-R) is a decreasing function is much wider on the critical line than in the 
chaotic phase; then, on the critical line the curves corresponding to a smaller A^ value are 



lower, while in the chaotic phase the contrary holds. As a consequence, if we average Y2{R) 
over R on the critical line, we get a quantity vanishing in the infinite size limit |10[, while 



the average weight of the attraction basins is finite and very close to the Random Map value 
in the chaotic phase 0. 

This difference and the non monotonic R behavior of 1^2 have a clear interpretation 
in the framework of the modular organization of Kauffman networks [Q. We thus postpone 
to that paper the discussion of our results. 
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Figure 1: The histogram of \ 2m — 1 |, where m is the average activity of a given node in 

a single cycle of a given network, is shown for K = 3 and N = 75. The peaks correspond 
to simple rational values of the argument: for instance, 1/8, 1/7, 1/6... 800 sample networks 
were generated, and 500 trajectories simulated on each of them. 
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Figure 2: Fraction s of nodes which are at time t in the same state in all the 200 different 
trajectories simulated. Data are average over 200 sample networks. K = 3, N = 50, 
100, 150 and 200 (N grows from top to bottom). The solid line is the mean field equation, 

s{t + l) = T.lo (f)^(t)^-'(l-^(t))' 
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Figure 3: Histogram of the networks with a fraction r of relevant elements. K = 3, N = 30 
(+), 40 (X), 50 and 60 (*). The networks generated were 20000, the relevant elements 
were identified simulating 200 trajectories on each of them. 
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Figure 4: The average period in networks with a fraction r of relevant nodes. K = 3, 
N — 30,40,50 and 60. Each curve was obtained generating at random 20000 networks and 
simulating 200 trajectories on each of them. 
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Figure 5: Probability density of the resettled number of unstttble elements Xu = U/N^^'^, where 
U is the number of unstable elements in a random network. The systems are at the critical 
point K — A, p — 1/A and system size ranges from 240 to 1600. 
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Figure 6: Probability density of the rescaled number of relevant elements Xr = R/N^^"^ , where 
R is the number of relevant elements in a random network. The systems are at the critical 
point K — A, p — 1/A and system size ranges from 120 to 1600. 




Figure 7: Average cycle length in networks with R relevant elements as a function of the 
resettled variable r — R/N^^^, for critical systems (K = A, p = 1/A) of size ranging from 120 
to 1600. 




Figure 8: Average value of I2 in networks with R relevant elements, as a function of the 
rescaled variable r — RsqrtN, for critical systems (K — 4, p = 1/4) of size from 120 to 
1600. 



